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ABSTRACT 

/>'  To  analyze  the  current  Aircraft  Battle  Damage  Repair  (ABDR) 
techniques,  a  finite  element  model  (FEM)  of  a  curved  composite  panel  with 
a  center  cutout  section  was  developed  and  verified  against  previously 
obtained  modal  test  results.  The  FEM  was  modified  using  different 
methods  to  exhibit  the  effects  of  riveting  on  an  aluminum  patch  to  repair 
the  section.  The  different  methods  were  contrasted  against  each  other  and 
previously  obtained  modal  test  results. 

The  panel,  a  '0/-45/+ 45/90],  layup  of  AS4/3501  graphite  epoxy 
measuring  12"  in  height  and  12.19"  across  its  arc  length  with  a  2"  by  4" 
cutout  at  a  0,  +45^  *45,  or  90  degree  orientation  was  modeled  using 
MSC/NASTRAN  QUAD4  elements.  The  first  repair  modeling  method  used 
a  mixture  of  the  aluminum  and  composite  panel  properties  to  represent  the 
overlap  region  of  the  patch.  The  second  method  used  ELAS2  and  CONM2 
elements  to  model  the  rivets,  and  modeled  the  panel  and  the  jpatch 
separately.  GAP  elements  prevented  the  patch  from  deforming  through  the 
panel. 

For  the  mixed  properties,  the  m  le  shapes  agreed  with  experimental 
data  at  lower  frequencies  but  were  distorted  at  the  higher  ones  and  were 
regularly  5-25%  high.  Modeling  the  rivets  individually  produced  matching 
mode  shapes  and  frequencies  regularly  within  5%  of  experimental  values.  ^ 


FINITE  ELEMENT  MODELING  OF  A  RIVETED  PATCH 
ON  A  COMPOSITE  PANEL 

I  INTRODUCTION 

The  many  advantages  of  the  material  properties  afforded  by  today’s 
advanced  composite  materials  have  brought  about  their  wide  use  in  the 
manufacture  of  many  of  today's  aerospace  vehicles.  But  the  application  of 
composites  to  existing  structural  problems  is  a  young  science,  and  not  a 
fully  explored  one.  Particularly  in  military  aircraft,  the  behavior  of 
damaged  composites  and  subsequent  repairs  are  of  great  interest. 

The  current  Aircraft  Battle  Damage  Repair  (ABDR)  techniques  for 
the  repair  of  battle  damaged  composite  panels,  based  on  simplicity,  speed, 
and  reliability  of  the  repair,  require  that  an  aluminum  patch  be  fastened 
with  bolts  over  the  damaged  section  of  the  composite  panel.  The  purpose 
of  ABDR  is  not  to  fully  restore  the  original  properties  of  the  panel.  It  is 
instead  intended  to  match  the  properties  closely  enough  that  the  aircraft  can 
be  returned  to  service  in  a  timely  manner  or,  if  damaged  beyond  the  repair 
abilities  of  the  unit,  allow  the  aircraft  to  sufely  fly  to  a  more  capable  repair 
facility  [18:1]. 


The  static  analysis  of  such  repairs  has  been  well  researched, 
including  the  development  of  the  McDonnel  Douglas  CREPAIR  program  for 
static  analysis  of  mechanically  fastened  repairs  to  composite  shells  [9:11]. 
However,  the  dynamic  analysis  of  a  mechanically  fastened  repair  is  a 
relatively  unexamined  though  important  problem. 

Understanding  of  how  these  repairs  might  affect  the  natural 
frequencies  of  the  panel  is  essential.  The  natural  frequency  is  the 
frequency  that  a  lightly  damped  structure  can  vibrate  at  a  continually 
increasing  amplitude  [20:85].  To  avoid  a  failure  of  the  section  under 
otherwise  normal  operating  condition,  the  effect  of  the  repair  on  the  natural 
frequency  should  be  ascertained.  This  provides  a  complicated  problem  for 
analysis:  determining  the  dynamic  behavior  for  a  curved  panel  with  a 
convoluted  geometry  and  complicated  material  properties.  Although 
closed-form  solutions  exist  for  many  .simpler  geometries,  the  combination 
of  unusual  geometric  boundary  conditions  and  complicated  material  analysis 
of  composite  shell  properties  present  great  difficiities  [15:205].  This 
problem,  however,  is  suitable  for  a  finite  element  oriented  solution. 
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finite  element  code  STAGS  C-1.  In  1986,  Cyr  [5]  followed  this  effort  by 
investigating  the  effects  of  orientation  of  a  2"  by  4*  interior  cutout  section 
on  the  modal  behavior  of  a  curved  rectangular  composite  panel  using 
interferometric  holography  techniques  and  the  finite  element  code  STAGS 
C-1.  In  1988,  Goodwin  [7]  analyzed  the  effect  of  riveting  an  aluminum 
patch  to  a  curved  rectangular  composite  panel  over  a  2"  by  4*  interior 
cutout  section  at  different  orientations  using  interferometric  holography 
techniques  and  the  finite  element  code  NISA.  Goodwin  was  limited  by  the 
capabilities  of  the  PC-based  NISA  code,  and  therefore  could  only  model  the 
symmetric  cases  and  modeled  the  repair  by  assuming  that  the  patch  was 
perfectly  bonded  to  the  panel.  The  perfect  bond  assumption  means  that  the 
panel  and  patch  maintain  the  identical  deformation  at  all  locations  where 
they  overlap.  In  1989,  Levraea  [12]  investigated  the  the  effects  of 
eccentrically  located  cutout  sections  on  the  modal  behavior  of  a  curved 
rectangular  composite  panel  subject  to  a  free-clamped-free-clamped 
boundary  condition  using  video  holography  and  STAGS-Cl. 

As  a  direct  percursor  to  the  attempt  to  model  the  individual  rivets  for 
this  thesis,  also  in  1988,  a  paper  by  Henriquez,  McLaughlin,  Montgomery, 
and  Ruggerio  [8]  discussed  the  study  of  modeling  of  various  thin  interfaces 
typically  encountered  in  small  resonant  structures,  including  two  surfaces 
bolted  together.  Two  free,  circular,  aluminum  plates  were  connected  with 
six  bolts  which  were  modeled  by  spring  elements  and  concentrated  mass 
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elements.  However,  the  a.ialytical  model  seemed  to  indicate  that  the  two 
plates  were  not  coupled  since  the  resonant  frequency  was  similar  to  that  of 
an  individual  plate  was  obtained  [8:430].  The  authors  did  not  attempt  to 
explain  why  this  might  have  occured. 

METHODOLOGY 

In  order  to  compare  to  Goodwin’s  experimental  results  and  to  avoid  a 
duplication  of  his  efforts,  the  curved  rectagular  composite  panel  used  in 
this  thesis  is  identical  to  the  one  described  in  Goodwin  [7:3. 1-3.1]. 
Similarly,  the  patch  and  the  fasteners  are  modeled  on  those  used  by 
Goodwin  [7:3.5,3.10]. 

Although  the  STAGS  C-1  analysis  was  based  on  the  Kirchhoff-Love 
hypothesis  and  neglected  t;ansverse  shear  while  the  NISA  analysis  included 
it  by  using  Reissner-Mindlin  theory,  both  codes  used  shell  elements  with 
higher-order  shape  functions  to  model  the  cuirved  panel  [6:2.22].  This 
investigation  used  MSC/NASTRAN  QUAD4  elements  that  combine  bending 
and  extensional  behavior  into  a  single  2-D  isoparametric  element  that 
includes  transverse  shear  effects.  The  curved  panel  was  then  modeled 
using  a  sufficent  number  of  flat  plate  element;  so  that  the  curvature  across 
each  individual  element  was  negligible.  This  allowed  the  use  of  simpler 
plate  elements,  as  well  as  a  simpler  plate  analysis  of  the  proi^erties  of  the 
composite  laminate.  The  QUAD4  elements  use  a  simpler  shape  function 


4 


than  the  NKTP  32  shell  elements  used  by  Goodwin,  and  provide  for  easier 
computational  solutions.  The  properties  of  a  composite  plate  are  dependent 
upon  the  curvature  of  the  mid-surface  of  the  plate,  and  a  composite  shell 
has  a  pre-existing  curvature  of  the  midsurface  prior  to  deformation  which  > 
greatly  increases  the  complexity  of  the  governing  equations. 

As  a  basis  for  comparison  to  Goodwin’s  results,  the  repair  was 
modeled  using  elements  with  a  blend  of  the  properties  of  the  composite 
panel  and  aluminum  patch,  including  the  ma^s  of  the  rivets,  in  the  region  of 
the  panel  where  the  patch  overlaps  the  panel  in  order  to  simulate  a  perfect 
bond  between  the  patch  and  the  panel.  The  main  effort  of  this  thesis  was  to 
more  accurately  model  the  effects  of  the  patch,  so  the  patch  and  panel  were 
modeled  as  separate  entities,  and  the  individual  rivets  were  individually 
modeled  to  connect  the  two  curved  surfaces. 

MSC/NASTRAN  ELAS2  spring  elements  were  used  to  model  the 
stiffness  contribution  of  the  rivets,  and  CONM2  concentrated  mass  elements 
were  used  to  represent  the  mass  of  each  rivet.  GAP  elements,  nonlinear 
elements  intended  for  use  with  nonlinear  geometries,  were  used  to  maintain 
the  proper  spacing  between  the  midsurfaces  of  the  panel  and  the  patch. 
MSC/NASTRAN  SINV  (Sturm  Inverse  Power  Method)  was  used  to  solve 
the  resulting  eigenvalue  problem. 

SDRC  I-DEAS  was  used  to  generate  the  geometry  of  the  models, 
define  the  physicaly  properties  of  the  elements,  write  the  MCS/NASTRAN 


data  deck,  and  to  graphically  represent  the  mode  shapes  resulting  from  the 
MSC/NASTRAN  ai.alysis. 

A  Sun4  workstation  was  used  to  run  the  programs  necessary  for  the 
execution  of  this  thesis.  SDRC  I-DEAS  Level  V  was  used  for  the  pre-  and 
post-processing,  and  MSC/NASTRAN  v66  was  used  to  perform  the 
analysis. 


II  THEORY 


OVERVIEW 

This  section  covers  the  three  theories  that  enabled  the  completion  of 
this  analysis:  Laminated  Plate  Theory  (LPT),  the  Finite  Element  Method 
(FEM)  and  the  Sturm  inverse  power  eigenvalue  solution  technique.  LPT 
allows  a.ply-by-ply  compilation  of  the  properties  of  a  laminated  composite 
to  develop  the  material  properties  of  the  laminate.  Using  FEM,  these 
material  properties  can  be  used  to  build  a  mathematical  model  of  the  curved 
panel.  The  dynamic  behavior  of  the  model  can  be  determined  by 
application  of  an  appropriate  solution  technique  to  the  eigenvalue  problem. 
For  this  thesis,  since  a  limited  number  of  eigenvalues  was  sought  within  a 
known  range,  the  Sturm  inverse  power  method  was  used. 


LAMINATED  PLATE  THEORY  (LPT) 


LPT  assume!,  for  simpUncalion  of  analysis,  that  a  perfect  bond 
exists  between  eacli  lamina,  that  each  lamina  can  be  considered  to  be  in  a 
ss,  and  that  displacement  is  continuous  in  the 
|.  LPT  combines  the  macromechanical  analysis  of 
with  these  three  assumptions  to  describe  the  layup  as  a 
whole.  This  analy  sis  is  briefly  covered  here  to  show  the  inclusion  of 
transverse  shear  stress  as  used  by  MSC/NASTRAN  [20:2.12-1  to  2.12-13]. 


state  of  plane  stres 
z-direction  [20:51] 
individual  laminae 


The  stress  is  related  to  strain  in  QUAD4  elements  by  the  relationship 

fj  ^4  ®  [e. 

<  in  '  «  G2  0  ■  2 

qj  0  0  t,Gj  [y 

where  {f}  =  [f,  f,  f^]T  the  membrane  forces/unit  leitgth; 

{m}  =  [m,  m,  m^]|  the  bending  moments/unit  length;  {q}  =  [q,  q^]!  the 
transverse  shear  force/unit  length;  {e}  =  [e,  €y  6,,]!  the  membrane  strains  at 
the  laminate  center;  {x}  =  [x»  X,  Xi,]T  the  curvatures;  and  {7}  =  [7,  7y][ 
the  transverse  shear  strains,  t,  is  the  effective  thickness  for  transverse  shear 


such  that  (t,  »  t‘5/6)  to  account  for  the  parabolic  distribution  of  shear 


stresses  through  a  laminate  of  thickness  t.  This  neglects  thermal  effects. 
For  a  laminate  with  N-piies  of  constant  thickness,  the  terms  [G|], 


[GJ  and  [G4]  are  defined  by 


M 

[0,1  -Z  /  [0.1, (b 

1  M- 
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H  ~  » 

[GJ  -  E  /  IG.],i*(b  (3) 
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where  the  [GJ  for  orthotropic  materials  is 


[GJ 


E, 


1-V,V,  l-VjV^ 
VjE,  Ej 


l-VjV,  l-V,Vj 


-  0 


0  a 


UJ 


(5) 


and  is  composed  of  the  material  properties  of  the  specific  lamina  in  the 
material  axis  system.  The  material  axis  system  is  transformed  to  the 
elemental  axes  by 

(OJ-[Tf[GJ[T] 


(6) 


where 


[T]« 


co66^  sm6^  cosBsinO 
sind^  co66^  >cos6sin6 
>2cos8sind  2cosdsm0  cos6^-sin9^ 


(7) 


The  [G3].  matrix  of  transverse  shear  moduli  is  deHned  by 
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[03].= 


On  Gi2 

I 

Ojj  Gjj 


which  is  transformed  to  the  element  coordinate  system  by 

[Gjl.  .  [WrtGjUWl 


(8) 


(9) 


where 


[W] 


cos6  sin6 
-sin6  cos6 


(10) 


The  value  of  the  transverse  shear  modulus  for  the  entire  laminate  is  derived 
from  the  transverse  shear  strain  energy  such  that  the  2  x  2  matrix  is 
defined  by 


N 


Ej,  k-l 


(11) 


where  En  and  E]}  are  the  1,1  and  2,2  terms  of  [GJ*  respectively  when 
[GJ*  is  the  same  as  [G2]  when  r  *■  0  and 


Rik  -  ♦  ;2k-i*k  ♦  i‘k*}*k*  02) 
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£imi£  ELEMENT  MEIim  (££M) 

FEM  takes  a  complicated  problem  and  divides  it  into  many  smaller 
problems  that  are  more  readily  solvable  [17:1].  Whereas  the  dynamic 
analysis  of  a  fully  clamped,  rectangular,  curved  composite  panel  is  difficult 
to  obtain  analytically,  the  solution  for  a  similar  simply  supported  flat  plate 
is  well  known  [22:87].  Similarly,  developing  the  stiffness  matrices  for  a 
curved  composite  shell  is  more  involved  than  the  flat  plate  condition  due  to 
the  initial  curvature  o^  the  midsurface.  However,  all  this  means  that  the 
solution  provided  by  the  FEM  is  only  an  approximation  of  the  exact 
solution,  and  is  not  exact.  It  can  be  shown  to  converge  on  the  exact 
solution,  given  enough  iterations  and  providing  that  compatibility 
relationships  are  satisfied  [3:126]. 

This  section  is  not  meant  as  a  thorough  development  of  Finite 
Element  Analysis  (FEA);  rather,  it  is  intended  to  highlight  those  portions  of 
the  theory  important  to  this  thesis.  First,  the  general  methods  and 
requirements  of  displacement  based  FEA  will  be  briefly  reviewed.  Then 
the  abilities  and  limitations  of  the  particular  MSC/NASTRAN  elements  used 
will  be  covered  in  the  section  most  applicable. 

FEA  TENETS.  The  first  step  in  the  creation  of  a  Hnite  element 
model  is  to  discretize  the  subject  [17:8].  For  this  project,  that  consisted  of 
dividing  the  curved  composite  panel  into  flat  elements  and  determining  the 
best  method  to  represent  the  effects  of  the  rivets.  Only  the  displacement 


(translational  and  rotational)  at  the  corners  or  ends  (nodes),  of  the  elements 
are  affected  by  the  displacement  functions.  They  are  related  to  each  other 
by  interpolation  formulae,  and  it  is  these  formulae  that  define  the  shape  and 
abilities  of  the  element.  These  shape  functions  are  arranged  into  a  matrix 
[N],  and  can  then  used  to  formulate  the  stiffness  matrix  for  the  element. 
First,  the  continuous  displacement  field  {u}  is  related  to  the  interpolation 
of  the  discretized  displacements  {d}  by 

(u)  -  [N]{d}  (13) 

Now,  the  discretized  displacement  can  be  related  to  the  elemental  strains  in 
much  the  same  manner  that  continuous  displacement  fields:  by 
differentiation.  The  equation 

(e)  »  [ajlu)  (M) 

for  a  set  of  continuous  displacements  yields 

(e)  -  (B]{d}  (15) 

for  a  set  of  discretized  displacements  where 

[B3  .  lailN]  (16) 

The  [d]  matrix  is  dependent  on  strain-displacement  relationships 
specific  to  the  problem  being  modeled. 
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As  can  be  found  in  almost  any  finite  element  text  (see  [3:110]  or 
[13:85-86]),  the  elemental  stiffness  matrix  [k]  can  be  derived  from  the 
principle  of  minimum  potential  energy,  resulting  in  the  volume  integral 

[k]  >  /l»f[E](B]dV  (17) 

V 

where  [E]  is  the  matrix  comprised  of  the  stress-strain  relationships  for  the 
discretized  system,  much  as  the  [B]  mainx  is  comprised  of  the  strain- 
displacement  relationships.  For  this  thesis,  [E]  is  the  stress-strain 
relationship  defined  in  Equation  1. 

It  is  important  to  note  that  the  [E]  and  [B]  matrices  determine  how 
fully  the  actual  physical  system  is  represented  in  the  finite  element  model. 
For  this  thesis,  a  full,  three-dimensional  model  was  not  created.  Rather, 
the  panel  was  modeled  using  idealized  two-dimensional  plate  elements 
which  are  considered  to  be  in  a  state  of  plane  stress.  Similarly,  the  full 
three-dimensional  effects  of  the  rivets  were  not  modeled;  they  wore 
represented  by  spring  elements  that  only  considered  one-dimensional 
effects.  These  simplifications  streamlined  the  modeling  process,  but  still 
resulted  in  an  accurate  solution. 

Since  every  differential  element  of  a  continuum  can  be  considered  to 
be  in  static  equilibrium  and  compatibility  prevails  everywhere,  a  anite 
element  model  should  attempt  to  reasonably  simulate  this.  As  stated  in 
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Cook,  the  extent  to  which  equilibrium  and  compatibility  conditions  can  be 
satisfied  are; 

1.  Equilibrium  of  nodal  forces  and  moments  is  satisfied. 

2.  Compatibility  prevails  at  nodes. 

3.  Equilibrium  is  usually  not  satisfled  across  interelement  boundaries. 

4.  Compatibility  may  or  may  not  be  satisfied  across  interelement 
boundaries. 

5.  Equilibrium  is  usually  not  satisfied  within  elements. 

6.  Compatibility  is  satisfied  within  elements. 

[124-126] 

Simply  stated,  these  mean  that  compatibility  and  equilibrium  must  be 
satisfied  between  nodes,  although  they  needn’t  necessarily  be  satisfied 
between  elements.  However,  the  incompatibilities  across  element  borders 
should  tend  toward  zero  as  more  e.lements  are  used.  This  is  important  when 
the  use  of  flat  elements  to  model  a  curved  panel  is  considered.  Modeling  a 
shell  with  plate  elements  provides  for  an  incompatibility  in  the  elemental 
out-of-plane  rotations,  as  is  depicted  in  Figure  1.  Because  the  elemental 


Figure  1. 

Z-axis  Rotation  Incompatibility 


z-axes  differ  by  the  small  angle  9,  neither  rotation  can  be  maintained  since 
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the  element  has  no  stiffness  in  that  direction  [14:181].  MSC/NASTRAN 
automatically  restrains  this  so  called  "drilling*  degree  of  freedom  in  order 
to  eliminate  this  incompatibility,  which  tends  to  artificially  stiffen  any 
model  using  elements  oriented  in  such  a  manner.  Increasing  the  number  of 
elements  should,  however,  decrease  this  effect. 

In  order  to  accomplish  a  dynamic  analysis  using  the  finite  element 
method,  not  only  must  the  stiffness  of  the  model  be  represented 
mathematically,  the  mass  of  the  model  also  must  be  expressed  in  matrix 
form.  Two  different  methods  are  available  to  accomplish  this  in 
MSC/NASTRAN,  the  lumped-mass  matrix  and  the  consistent  mass  matrix. 

The  lumped  mass  method  simply  apportions  the  mass  of  the  element 
equally  among  the  translational  degrees  of  freedom.  This  assumes  that  the 
material  within  the  mean  location  on  either  .^ide  of  the  particular 
displacement  behaves  like  a  rigid  body  while  the  rest  of  the  element  does 
not  participate  in  the  motion.  Because  this  excludes  the  dynamic  coupling 
that  exists  between  the  elemwnt  displacements,  it  has  the  advantage  of 
creating  a  diagonal  matrix;  but,  it  usually  ignores  the  effects  of  rotary 
inertia  [17:365].  For  most  dynamic  cases  where  the  displacements  are 
small,  the  effects  of  rotary  inertia  are  not  notable,  so  neglecting  them 
should  not  be  a  cause  for  error. 

The  consistent  mass  matrix  is  defined  by  the  volume  integral 


(18) 


[ml  p  / [NfpCNldV 

V 

where  p  is  the  density  of  the  element  and  [N]  is  the  shape  function.  Since 
the  consistent  mass  matrix  depends  on  the  shape  function,  it  will  generally 
be  a  fully  populated  matrix  [13:475].  The  consistent  mass  matrix 
automatically  includes  the  rotary  inertia  terms,  and  the  off-diagonal  terms 
represent  the  dynamic  coupling  terms.  This  method  involves  an  integral, 
which  must  be  approximated  for  use  on  a  computer,  and  is  therefore  more 
numerically  complex  than  the  lumped  mass  matrix  approach,  and  requires 
more  computer  storage  space.  Also,  the  diagonal  lumped  mass  matrix  can 
facilitate  the  computational  procedures  [17:365].  For  this  thesis,  though, 
the  results  differed  by  less  than  1%  depending  on  the  mass  matrix  used. 
Since  the  lumped  mass  matrix  provided  for  a  faster  solution,  it  was  applied 
to  all  the  model  solutions  for  this  project. 

QUAD4  ELEMENT.  The  MSC/NASTRAN  QUAD4  e’ement  was 
developed  using  a  four-noded  isoparametric  formulation  so  that  it  could  be 
used  for  an  arbitrary  quadrilateral  geometry.  Its  development  is  well 
documented  and  only  the  highlights  are  presented  here  [14:175].  The 
element  was  developed  for  use  with  anisotropic  properties,  and  features  an 
elastic  coupling  between  bending  and  stretching  and  an  augmentation  of 
transverse  shear  flexibility.  To  do  this  requires  modifying  the  standaid 
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isoparametric  theory  for  both  membrane  and  bending  action. 

The  Hrst  change  involves  computing  the  shear  strain  at  the  center 
of  the  element  (^  -  v  ~  0).  This  value  is  used  to  evaluate  stresses  and 
strain  energy  at  a  2  x  2  array  of  Gauss  points,  together  with  strains  6.  and 
computed  directly  at  these  points  [14:177].  This  is  a  modiflcation  to  the 
standard  isoparametric  theory  for  membrane  action,  and  permits  element 
edges  to  be  curved. 

To  modify  the  bending  action,  the  out-of-plane  displacements  are 
now  defined  as  a  double  power  series  in  (  and  ij: 

W  ♦  .  .  .  (19) 

If  the  coefficients  are  arranged  in  a  pyramid  as  shown  in  Figure  2, 


Figure! 

Coefficient  Triangle 


then  the  first  two  row.  represent  the  rigid  body  motions  and  constant 
transverse  shear  strains,  the  third  row  represents  the  constant  curvatures, 
and  the  fourth  row  represents  curvatures  that  vary  linearly  with  (  and  if. 
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The  W20  and  Wq,  terms  require  reduced  order  integration  of  the  energy  due 
to  the  transverse  shear,  and  the  w,o  and  Wo3  terms  can  be  corrected  by  the 
introduction  of  residual  bending  flexibility.  This  is  done  by  setting  the 
conditions  for  curvature  dependence  to  be  the  curvature  compatibility 
conditions  for  a  plate  with  zero  transverse  shear  strains. 


.  2^ 

dx 

Sy 

.  2^ 

dx 

which  are  satisfied  if  the  twist  used  compute  the  strain  energy  is 

where  x!y  is  the  value  of  twist  computed  at  any  integration  point  g  and  xi, 
is  the  value  of  twist  computed  at  {  «  i;  •  0.  Thus  the  energy  due  to 
bending  is  computed  by  a  numerical  integration  at  a  2  x  2  array  of  Gauss 
points.  However,  this  leads  to  an  energy  deficiency.  It  is  correctable, 
though,  by  the  introduction  of  residual  bending  flexibility. 

ELAS2  ELEMENT.  The  ELAS2  element  (shown  in  Figure  3)  is 
simply  a  three-dimensional  translation  spring  element.  It  can  act  as  up  to 
three  independent  translational  springs  acting  along  three  orthogonal 
directions,  with  each  direction  having  its  own  spring  stiffness  constant.  Its 
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Figure  3. 

ELAS2  Element 

stiffness  matrix  is  simply  defined  by 

f^  -  (22) 

where  Ki  is  the  directional  stiffness,  Xi  is  the  directional  displacement,  and 
fi  is  the  force  acting  on  it.  It  was  used  in  this  project  to  simulate  the 
stiffness  of  the  rivets  by  equating  each  Ki  with  Young's  Modulus  of 
Elasticity  of  aluminum. 

GAP  ELEMENT.  The  GAP  element  (shown  in  Figure  4)  has  a 


coaiditional  stiffness  matrix.  With  the  definition  of  an  initial  gap  opening, 
the  GAP  element  will  exhibit  different  values  of  axial  stiffness,  depending 
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on  whether  the  initial  separation  distance  has  increased  or  decreased.  It 
can  also  provide  or  transverse  stiffness  when  the  gap  separation  is  open 
(separation  greater  than  initially  defined),  as  well  as  coefficients  of  friction 
in  transvers  e  directions.  These  latter  capabilities  were  not  explored  in  this 
thesis,  as  nv  transverse  motion  was  expected.  The  GAP  element  was  used 
to  maintain  the  initial  separation  between  the  midsurfaces  of  the  panel  and 
the  patch.  This  will  be  more  fully  explained  in  the  next  chapter. 

EIGENYALUE  AMALXSIS. 

From  Newton’s  Second  Law,  the  following  equation  governs  motion 
for  a  discretized  system  [21:144]: 

[M]{1)  ♦  [C](i}  ♦  [K){x}  -  {F(t)}  (23) 

where  [M]  is  the  mass  matrix,  [C]  is  the  damping  matrix,  [K]  is  the 
stiffness  matrix,  {x}  is  the  position,  and  {F(t)}  is  the  external  forces  acting 
on  the  system.  For  the  free  vibration  of  a  thin  plate  where  damping  caii  be 
considered  negligible,  [C]  -  0  and  {F(t)}  *  0.  Now  the  governing 
equation  is 

(Mid)  ♦  [Klix)  -  (0)  (24) 

By  assuming  a  displacement  function  of  the  form 
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(x)  ■ 


(25) 


and  substituting  into  Equation  (24),  the  characteristic  eigenvalue  problem  is 
derived: 

[IKI  -  -  0  <26) 

for  which  non-trivial  solutions  exist  only  when 

|[K]  -  «^M]|  -  0  <27) 

The  Sturm  Inverse  Power  Method  (SINV),  which  is  a  variant  of  the 
inverse  power  method,  was  used  to  solve  for  the  eigenvalue  \  =  using 
shifts  for  faster  convergence,  sweeping  to  eliminate  effects  of  prior 
eigenvalues,  and  the  application  of  Sturm  sequence  polynomials  to  aid  the 
control  of  the  iteration  process  [20:4.4-40]. 

Shifting  the  solution  algorithm  to  an  initial  A,  such  that 

(28) 

so  A  replaces  X  as  the  eigenvalue  parameter,  producing 

(K-  A,M](*)  -  A[Ml{(i)}  (29) 

which  is  the  basis  of  the  iteration  scheme 
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[K-A^M]{4>.}  -  [M]{v,} 


(30) 


{u.}.±{4)^}  (31) 

[20:4.4-24] 

where  C,  is  equal  to  the  element  of  with  the  greatest  absolute  value  and 
1/Ca  converges  to  A|,  the  shifted  eigenvalue  nearest  the  shift  point,  and 
{u.}  converges  to  the  corresponding  eigenvector. 

This  eigenvector  can  be  eliminated  via  a  sweeping  procedure,  since 
all  the  eigenvectors  are  orthogonal  [15:119].  For  a  discretized  system,  the 
number  of  degrees  of  freedom  equals  the  number  of  eigenvalue  solutions  to 
Equation  (27)  [17:153].  Therefore,  once  the  Hrst  eigenvalue  has  been 
found,  the  size  of  characteristic  equation  can  be  reduced  to  accelerated  to 
determination  of  the  next  eigenvalue. 

To  solve  Equation  (30),  the  SINV  solution  method  decomposes  the 
[K  -  X,M]  matrix  into  triangular  factors  such  that 

(K-X.Ml  .  (LHDJtLf  (32) 

where  [L]  is  a  lower  triangular  matrix  with  ones  on  its  diagonal  aiid  [D]  is 
a  diagonal  matrix  consisting  of  the  terms  on  the  diagonal  of  the  [K  -  X^M] 
matrix.  Then  the  characteristic  polynomials  of  the  leading  principal 
submsltrix  of  the  [K  •  X,M]  matrix  can  be  written  as  a  Sturm  sequence  of 
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polynomials.  It  can  then  bw  shown  that  the  number  of  sign  agreements  in 
the  Sturm  sequence  corresponds  to  the  number  of  negative  values  in  the 
diagonal  [D]  matrix,  which  equals  the  number  of  eigenvalues  of  [K  -  X^M] 

below  X,  [20:4.4-34].  This  is  calculated  for  each  shift,  and  is  used  in  the 
control  of  the  iteration  process. 


23 


Ill  FINITE  ELEMENT  MODELING 


In  order  to  validate  this  analysis,  three  steps  were  necessary.  First 
the  properties  of  the  [0/-45/+ 45/90],  AS4/3S01  Hercules  graphite  epoxy 
needed  to  be  modeled.  Then  a  model  of  the  full  undamaged  panel  needed 
to  be  created  and  solved,  and  the  results  compared  to  the  experimental  data 
of  Goodwin.  Finally,  four  more  models  were  necessary,  one  for  each  of 
the  following  cutout  section  orientations:  Ot  90?  +45?  and:  -45?  The 
results  from  these  analyses  were  then  compared  to  the  experimental  results 
of  Goodwin. 

Once  the  validation  of  the  panel  modeling  was  complete,  the  cutout 
section  models  were  changed  by  adding  an  aluminum  patch  over  the  cutout 
section.  This  was  done  by  two  methods.  The  first  assumed  that  the  patch 
was  peifectly  bonded  to  the  panel.  The  other  method  modeled  the  rivets 
individually. 

LAMINATE  ERQPERIX  MOPEUNG 

The  orthotropic  material  properties  used  by  Goodwin,  listed  in 
Table  1 ,  were  used  to  create  an  orthotropic  material  table  in  SDRC 
I-DEAS.  Using  the  Laminate  Modeling  Task  to  define  the  orientation  and 
thickness  of  each  ply  allowed  SDRC  I-DEAS  to  calculate  the  laminate 
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Table  1.  AS4/3S01  Material  Properties 


1  Material  Property 

Value  1 

Modulus 

of 

Elasticity 

E, 

18.84e6  psi  | 

E, 

1.468e6  psi  | 

Shear 

Modulus 

Gi2 

0.9099e6psi  | 

0-23 

0.9099e6  psi  | 

G„ 

0.7279e6  psi  | 

Poisson’s 

Ratio 

*'12 

0.2799 

•'ll 

0.0218  , 

Density 

fi 

1.516e-4* 

Thickness 

t 

0.005"  1 

(*  units  for  p  are  lbf*sVin^) 


properties  which  are  displayed  graphically  in  Figure  5.  These  material 
properties  were  then  oriented  on  all  panel  elements  so  that  the  primary 


element  material  axis  was  parallel  to  the  global  z-axis. 
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Figure  5.  Laminate  Properties 
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FULL  PANEL  MODEL 


The  full  panel  was  modeled  first  by  using  elements  approximately 
Vi"  square  (0.500*  x  0.508")  as  is  shown  in  Figure  6.  The  outer 


Figure  6.  Full  Panel  Model 


dimensions  (12"  high  by  12.19"  arc  length  across  a  12"  radius  of  curvature) 
are  the  same  for  all  the  models,  and  the  same  as  those  used  in  experiments 
of  Goodwin,  Cyr,  and  Walley. 

To  insure  that  the  model  would  convergence,  the  mesh  was  refined  so 
that  the  elements  were  approximately  14”  square  (0.250”  x  0.254”)  as  is 
shown  in  Figure  7.  These  results  were  contrasted  against  the  frequency 


Figure  7.  ReHned  Panel  Model 

results  for  the  first  five  modes  of  the  previous  model  and  the  experimental 
data,  as  well  as  the  cpu  time  required  to  compute  the  solution  in  order  to 
determine  the  optimal  modeling  option  (see  Table  2).  On  this  basis,  it  was 

Table  2.  Refinement  Results 


Source 

ft 

cpu  (s) 

V4"  Mesh 

558 

591 

761 

825 

916 

3359.3 

14"  Mesh 

550 

571 

755 

820 

905 

15918.3 

Goodwin  Exper. 

514 

519 

716 

727 

821 

N/A 

(*  All  frequencies  in  hert2) 

decided  to  use  the  Hrst  (coarser)  mesh  as  a  basis  for  construction  of  the 
other  models,  as  the  increase  in  accuracy  was  severely  outweighed  by  the 
increase  in  computational  effort. 


cmmx 


Q°AND  90°  SECTIONS.  These  models,  shown  in  Figure  8  and 
Figure  9,  were  very  similar  to  the  fv’i  panel  model  since  the  c.^tout  sections 


needed  to  be  refined  within  0.75"  of  the  cutout,  to  a  0.25"  x  0.25"  grid. 
Triangular  7.RIA3  elements  were  used  to  transition  from  the  outer  to  the 
inner  grid.  Although  MSC/NASTRAN  advocates  the  use  of  RSPLINE 
elements  for  this  task,  doing  so  would  disallow  the  use  of  SDRC  I*DEAS  as 
a  post-processor,  since  that  element  is  not  supported  by  I-DEAS  Level  V. 

*45.“  SECTIONS.  Since  the  mode  shapes  of  a  cylindrical  section 
tend  to  be  aligned  with  the  global  z-axis  of  the  panel,  these  models 
presented  a  set  of  contradictory  demands;  the  elements  should  be  aligned 
with  cutout  section  for  the  most  physically  accurate  model,  yet  the  elements 
should  be  aligned  with  the  expected  mode  shapes  so  that  the  interpolation 
functions  are  best  able  to  approximate  the  actual  displacement  function. 
Initial  attempts  to  align  the  elements  with  the  cutout  were  unsuccessful,  as 
the  mode  shapes  were  deformed  and  in  some  cases  unrecognizable.  Using  a 
34  X  34  alleviated  this  problem,  though.  Its  elements  were 
0.353"  X  0.359"  with  .a  diagonal  measuring  0.503?  Thus,  a  cutout  section 
measuring  3.95*  x  1.97*  was  created,  and  the  elements  were  aligned  with 
the  expected  displacement  functions.  These  models  are  shown  in  Figure  10 
and  Figure  11. 
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Figure  10.  +45 "  Cutout  Section  Model 


Figure  11.  -45*  Cutout  Section  Model 
The  presence  of  triangular  elements  along  the  border  of  the  cutout 
sections  should  not  present  any  difFiculties  in  achieving  an  accurate 
solution,  since  the  cutout  section  will  not  create  the  type  and  magnitude  of 
stress  concentrations  in  a  dynamic  analysis  as  it  would  in  a  static  loading 
analysis. 
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PAICHMOPEUNq 

PERFECT  BOND  MODEL.  Assuming  that  the  patch  is  perfectly 
bonded  to  the  panel  has  been  the  standard  assumption  in  prior  attempts  to 
model  the  dynamic  behavior  of  a  repaired  panel.  This  was  done  in  this 
project  to  provide  a  baseline  for  comparison  to  the  individual  rivet 
modeling  method.  To  prepare  this  model,  the  mesh  was  divided  into  three 
regions:  the  panel,  the  patch,  and  the  overlap  (see  Figure  12).  The  panel’s 
characteristics  were  unchanged  from  the  cutout  model,  and  the  patch  was 
modeled  using  the  isotropic  aluminum  properties  shown  in  Table  3.  The 


Table  3.  Patch  Material  Properties 


I  Material  Property 

Value 

Modulus  of  Elasticity 

E 

10.5e6  psi 

Poisson’s  Ratio 

9 

0.33 

Density 

P 

2.6i7e*4' 

Thickness 

t 

0.032- 

(*  units  for  p  are  lb,*s*/in*) 


overlap  region  was  modeled  using  a  mixture  of  the  aluminum  and  laminate 
properties,  which  were  calculated  using  a  method  similar  to  the  rule  of 
mixtures  used  to  determine  the  properties  of  a  composite  based  on  the 
flber/matrix  volume  ratio.  In  order  to  be  able  to  do  this,  equivalent 
transverse  isotropic  properties  were  needed  for  the  laminate.  Since  the 
layup  is  balanced  and  symmetric,  the  following  equations  apply: 
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Table  4.  Equivalent  Laminate  Material 
Properties 


1  Material  Property 

Value 

Modulus 

of 

Elasticity 

E. 

7.515e6  psi 

E, 

7.515e6  psi 

Poisson’s 

Ratio 

0.293 

0.293 

Shear  Modulus 

2.906e6  psi 

The  resulting  material  properties  for  the  overlap  section  can  be  found 


in  Table  5. 


Table  5.  Overlap  Material  Properties 


1  Perfectly  Bonded  Patch/Panel  Material  | 

1  Material  Property 

Value 

Modulus  of  Elasticity 

E 

8.8e6  psi 

Poisson’s  Ratio 

p 

0.309 

Density 

p 

2.01 8e-4’ 

(*  units  for  p  are  Ibf’sVin*) 


&1V 


EX  MODELING.  The  ABDR  requirements  dictated  the  size  of 


patch  and  I 
generated 
panel,  the 


the  placement  of  the  rivets  (see  Figure  12).  The  patch  whs 
as  a  reHned  square  3.3*  x  5.3*  mesh  0.036*  above  the 
distance  from  the  midsurface  of  the  panel  to  the  midsurface  of 


the  patch. 


Spring  elements,  lumped  mass  elements,  and  gap  elements  were 


Figure  12.  Rivet  Placement 


used  to  model  the  rivets  and  separation,  as  can  be  seen  in  Figure  13. 
Node-to-node  translational  spring  elements  were  generated  at  the  rivet 
locations  in  SDRC  I-DEAS,  which  translated  into  ELAS2  elemenU  in 
MSC/NASTRAN.  SDRC  I-DEAS  lumped  mass  elements,  which  became 
CONM2  MSC/NASTRAN  elements,  were  generated  at  both  end  nodes,  each 
representing  half  of  the  rivet  mass  (2.139e-6  Ibf’sVin).  SDRC  I-DEAS  gap 
elements,  equivalent  to  MSC/NASTRAN  GAP  elements,  were  applied  to  the 
nodes  on  the  edge  of  the  cutout  section  of  the  panel  and  the  outer  edge  of 
the  patch,  and  the  opposing  nodes  of  the  opposite  surface.  The  nodes 
between  the  edges  did  not  have  gap  elements  separating  the  two  surfaces 
since  this  was  found  to  over-constrain  the  model  and  prevent  solution. 
However,  the  separation  of  these  inner  nodes  did  not  vary  in  any 
appreciable  manner,  due  largely  to  the  inherent  stiffness  of  the  elements. 
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Figure  13.  Close  View  of  Patch 


In  Figure  13,  the  bottom  layer  of  quadrilateral  elements  is  the 
composite  panel,  and  the  top  layer  of  quadrilateral  elements  is  the 
aluminum  patch.  The  lumped  mass  elements  are  represented  by  squares 
around  their  nodal  locations,  and  the  spring  elements  show  the  rivet 
placement  and  modeling.  Along  the  outer  edge  of  the  patch  are  the  GAP 
elements,  as  well  as  along  the  edge  of  the  cutout  (in  the  top  left  corner). 
The  use  of  GAP  elements  makes  this  model  work.  An  initial  solution  was 
run,  without  the  presence  of  GAP  elements.  The  patch  and  panel  vibrated 
as  separate  entities,  connected  only  at  the  rivet  locations,  and  often 
deformed  in  opposite  directions,  passing  through  one  another.  This  may  be 
what  happened  to  the  model  described  by  Henriquez,  et  al  [8].  Since  their 
paper  did  not  fully  describe  the  model  or  the  results  or  even  show  the  mode 
shapes,  but  did  indicate  that  the  two  bolted  plates  did  not  appear  to  be 
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coupled,  this  seems  likely. 

SDRC  I-DEAS  was  used  to  produce  the  MSC/NASTRAN  bulk  data 
deck  and  case  control  deck.  Unfortunately,  I-DEAS  Level  V  cannot 
consistently  label  continuation  notation  for  CELAS2  cards,  and  does  not 
write  the  orientation  vector  properly  to  the  CGAP  cards.  This  required 
manually  verifying  the  continuation  notation  for  all  the  CGAP  cards,  and 
correcting  them  when  necessary.  To  fix  the  CGAP  cards,  a  unit  vector  was 
defined  in  the  radial  direction  of  the  cylindrical  coordinate  system  to 
provide  orientation  for  the  GAP  elements  and  added  to  each  CGAP  card. 


IV  BS^UXIS 


This  chapter  will  summarize  and  discuss  the  resulting  solutions  of  the 
finite  element  models  in  three  sections.  First,  the  full  panel  will  be 
compared  to  the  experimental  results  of  Goodwin  and  contrasted  against 
prior  finite  element  analysis  using  NISA.  Next,  the  cutout  section  models 
will  be  compared  to  the  experimental  results  of  Goodwin  and  contrasted 
against  prior  Hnite  element  analysis  using  NISA.  Finally,  the  individual 
rivet  model  will  be  compared  to  Goodwin’s  experimental  results  and 
contrasted  against  the  perfect  bond  model  and  the  prior  NISA  results.  Note 
that  a  discussion  of  the  effects  of  the  cutout  sections  on  the  mode  shapes  of 
the  panels  is  not  included,  since  both  Goodwin  [7]  and  Cyr  [4]  provided 
detailed  explanations. 

The  mode  shapes  are  shown  for  each  model  using  displacement 
contours.  Although  SDRC  I-DEAS  is  capable  of  producing  deformed 
geometries  from  the  MSC/NASTRAN  output  Hies,  displacement  contours 
are  closely  related  to  the  holographic  mode  shapes  provided  by 
experimental  analysis  snd  therefore  provide  a  better  basis  for  comparison. 

FULL  EANEL  MQD£].. 

This  model,  as  previously  stated,  was  created  to  insure  that  the 


modeling  procedures  used  in  this  analysis  represent  an  accurate 
approximation  of  the  actual  model.  The  frequency  results,  which  can  be 
seen  in  Table  6,  show  that  this  model  is  artificially  stiff.  This,  as  was 


Table  6.  Full  Panel  Results 


_ m 

U 

r. 

U 

f, 

Full  Panel  Model  (PPM) 

55S 

591 

761 

825 

916 

Ooodwin  Experinentel  (OE) 

514 

519 

716 

727 

821 

%  DifTerenee  PPM  vt  GE 

•»-S.56S 

•flS.t* 

+6.2t% 

•t-13.5S 

•l•11.6« 

Befined  PPM 

550 

571 

755 

820 

905 

%  Difference  vt  GE 

+7.00* 

-flO.OX 

-t-5.45S 

4-12.8K 

+  10.2* 

NISA  (Goodwin) 

561 

5tt 

797 

846 

N/A 

%  Difference  vt  GE 

4-13.3S 

•fll.3S 

■t-16.4« 

C  aU  fr«4U«aeiM  ia  btitz) 


Stated  in  the  theory  section,  is  due  to  the  automatic  restraining  of  the 
elemental  z-axis  rotation  to  insure  compatibility.  It  is  evident  that  refining 
the  model  improved  the  accuracy  of  the  analysis,  but  doing  so  increased  the 
cpu  time  required  for  the  solution  by  nearly  five  times  (3359.3  vs  15918.3 
cpu  seconds).  The  mode  shapes  for  both  models  are  nearly  identical,  as 
can  be  seen  in  Figures  14  through  23.  The  mode  shapes  also  match  those 
reported  by  Goodwin  in  both  his  experimental  holography  work  and  NISA 
finite  element  wdrk.  However,  both  the  original  and  the  reHned  full  panel 
inodels  show  closer  agreement  with  the  experimental  natural  frequencies 
than  the  NISA  results. 
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Figure  18.  Full  Panel  Model 
Mode  3,  f  =  761  hz 


Figure  19.  Refined  Full  Panel  Model 
Mode  3,  f  »  755  hz 


CUTOUT 


MODEL  RESULTS 


CUTOUT  SECTION  MODEL.  The  natural  frequencies  of  this 
model  can  be  found  in  Table  7.  They  are,  as  expected,  higher  than 

Table  7.  0°  Cutout  Section  Model  Results 


fT 

U 

f, 

L 

u  1 

1  0*  Cutout  Model  | 

491 

560 

662 

758 

812  1 

1  Goodwin  Exper.  (GE)  | 

445 

507 

606 

708 

749 

1  %  Difference  | 

+  10.:* 

+  10.4* 

+9.24* 

+  7.06* 

+  8.41* 

1  NISA  (Goodwin)  | 

499 

581 

681 

789 

N/A  1 

1  51  Diff.  vt  GE  1 

+  12.1* 

+  14.6* 

+  12.4* 

+  11.4* 

...  1 

(*  all  frequencies  are  is  hertz) 


Goodwin’s  experimental  results  and  more  accurate  than  Goodwin’s  NISA 
results.  The  mode  shapes  are  shown  in  Figures  24  through  28.  Mode 
switching  is  experienced  with  the  second  and  third  mode  shapes,  relative  to 
Goodwin’s  experimental  work.  However;  these  results  agree  with 
Goodwin’s  NISA  work,  as  well  as  Cyr’s  experimental  work.  As  Goodwin 
noted,  mode  switching  is  a  phenomenon  wherein  a  symmetric  and  an 
adjoining  asymmetric  mode  shape  exchange  frequencies.  It  is  felt  in  this 
instance,  though,  that  Goodwin  experienced  mode  switching  in  his 
experimental  work,  and  that  the  order  of  mode  shapes  presented  here  is  the 
proper  one.  Also,  the  second  and  third  mode  shapes  of  this  model 
correspond  to  the  same  respective  mode  shapes  of  the  full  panel  model. 


Figure  29.  0*  Cutout  Effects 


experimental  work  of  Goodwin,  and  the  improvement  over  the  NISA  results 
is  also  obvious. 

CI/TQUT  SECTION  MQPELS 

CUTOUT  SECTION  MODEL.  The  results  of  this  model  do  not 
show  perfect  agreement  with  Goodwin’s  experimental  results.  The 
frequency  results  show  that  once  again  the  natural  frequencies  of  this  model 
follow  the  trend  of  the  natural  frequencies  of  the  experimental  work,  only 
higher.  The  mode  shapes  only  match  Goodwin's  results  for  the  first  four 
natural  frequencies,  although  mode  switching  occurs  between  the  first  and 
second  mode  shapes.  However,  this  model's  first  two  mode  shapes  agree 
with  those  from  Goodwin's  NISA  work,  as  well  as  the  first  two  mode 
shapes  resulting  from  Cyr's  experimental  and  STAGS  C>1  work.  Therefore 
the  first  four  mode  shapes  were  immediately  accepted  as  accurate. 

The  sixth  mode  shape  of  this  model,  though,  matches  Goodwin’s  fifth 
mode  shape  while  the  fifth  mode  matches  nothing.  Since  Goodwin’s  fifth 
mode  matches  his  NISA  work,  as  well  as  both  the  experimental  and 
STAGS  C*1  FEA  work  of  Cyr,  the  initial  reaction  was  to  regard  this  new 
mode  shape  as  a  repercussion  of  this  model's  use  of  plate  elements  to  model 
the  curved  panel.  The  mesh  was  reHned  and  the  solution  rerun  in  order  to 
ascertain  if  increasing  the  number  of  degrees  of  freedom  of  the  model 
would  resolve  the  situation.  The  reHned  mesh  frequencies  results  were 


indeed  more  accurate,  as  can  be  seen  in  Table  8,  but  the  mode  shapes  did 


Table  8.  90"  Cutout  Section  Model 


fT 

U 

f4 

fj 

90*  Cutout  Model 

578 

582 

716 

765 

,  779 

Goodwin  Exp.  (GE) 

519 

546 

656 

712 

733 

%  Difference 

-l-11.4ft 

•»-6.60ft 

+9.15ft 

+  7.44  ft 

+6.28ft 

Refined  Mesh 

570 

575 

709 

744 

782 

%  Diff.  vs  GE 

+9.83ft 

+5.31ft 

+  8.08ft 

+5.20ft 

+  3.96* 

NISA  (Goodwin) 

584 

585 

795 

812 

N/A 

ft  Diff.  vs  GE 

-t-12.5ft 

+7.16ft 

+  21. 2ft 

+  14. Oft 

— 

(*  all  frequencies  are  in  hertz) 


not  change  appreciably,  as  can  be  seen  in  Figures  30  through  41. 

As  a  consequence  of  the  results  of  the  0"  repaired  models  and  the 
±45"  cutout  section  models  (where  similar  situations  occurred),  this 
behavior  is  now  thought  to  be  merely  mode  switching  between  this  model's 
fifth  and  sixth  modes.  In  experimental  analysis,  mode  switching  may  be  a 
function  of  the  slightest  difference  in  testing  conditions,  whether  it  is  a 
slightly  different  boundary  condition  or  a  material  imperfection  in  the  test 
object  itself  [7:6-11].  For  the  finite  element  model,  the  higher  natural 
frequencies  are  more  sensitive  to  the  boundary  conditions  [19:232-233]. 
The  artificial  restraining  of  the  elemental  z-axes  for  this  model  could 
therefore  be  responsible  for  the  mode  switching  at  higher  frequencies. 
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Figure  30.  90**  Cutout  Section  Model 
Mode  1,  f  *  578  hz 


Figure  31.  90*  Cutout  Section  Refined  Model 
Mode  1,  f  »  570  hz 


Figure  32.  90**  Cutout  Section  Model 
Mode  2,  f  =  582  hz 


Figure  34.  90°  Cutout  Section  Model 
Mode  3,  f  =  716  hz 


Figure  35.  90°  Cutout  Section  ReHned  Model 
Mode  3,  f  »  709  hz 
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Figure  40.  90"*  Cutout  Section  Model 
Mode  6,  f  =  802  hz 


Figure  41.  90**  Cutout  Section  Refined  Model 
Mode  6,  f  a*  790  hz 


For  this  cutout  orientation,  the  first  natural  frequency  is  higher  than 
that  of  the  full  panel.  Goodwin  also  noted  this  behavior,  although  he  noted 
it  the  first  two  natural  frequencies.  As  shown  graphically  in  Figure  42, 


Figure  42.  90**  Cutout  Effects 


this  model  exhibited  nearly  the  same  frequency  at  the  second  mode  for  both 
the  90°  cutout  section  and  the  full  panel  model.  Since  the  natural 
frequencies  are  directly  proportional  to  stiffness  and  inversely  proportional 
to  mass,  the  behavior  of  this  model  seems  to  indicate  that  more  mass  was 
lost  than  stiffness  with  the  90°  cutout  section  [6:3032]  for  the  fii'st  mode 
shape.  This  makes  intuitive  sense,  since  the  greatest  stiffness  is  provided 
by  the  circumferential  fibers  (the  90°  fibers)  and  the  fewest  would  be 
damaged  by  this  cutout.  Furthermore,  the  mode  shapes  become 
successively  more  complex,  which  would  emphasize  the  effects  of  the  lost 
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stiffness.  Because  of  the  mode  switching,  Goodwin’s  second  mode  shape 
was  simpler  than  his  first.  This  would  explain  the  respectively  lower  result 
for  the  second  natural  frequency  of  this  model. 

+41“  CUTOUT  SECTION  MODEL.  This  model  displays  the  same 
characteristics  as  the  90“  cutout  section  model,  wherein  the  first  four  mode 
shapes  match  Goodwin’s  results  but  the  sixth  mode  shape  of  this  model 
matches  Goodwin’s  fifth.  Similarly,  this  model’s  Hfth  mode  matches  that 
found  by  Cyr  in  both  his  experimental  and  STAGS  C-1  FEA  work.  The 
natural  frequencies  of  this  model,  shown  in  Table  9,  are  higher  than 


Table  9.  +45“  Cutout  Section  Model  Results 

I  r,  f,  I  f,  f. 


+  45*  Cutout  Model 

1 

592 

711 

753 

823 

Goodwin  Exper. 

1  soo 

538 

631 

685 

775 

%  Difference 

1  +11.2X 

+  10.0* 

+  12.7* 

+  9.93* 

+  6.19* 

(*  all  frequencies  ere  in  hertz) 

Goodwin’s  experimental  results,  as  expected.  The  mode  shapes  are  shown 
in  Figures  through  48.  The  discrepancy  between  this  model’s  fifth  and 
sixth  modes  is  considered  to  be  modes  switching,  especially  in  light  of  the 
results  of  the  -45“  cutout  section  model.  The  sixth  mode  of  the  -45“  cutout 
section  model  is  the  transpose  of  the  +45*  cutout  section  model’s  fifth 
mode  shape  and  matches  the  fifth  mode  shape  predicted  by  Cyr’s 
STAGS  C-1  work.  This  will  be  more  fully  explained  in  the  next  section. 


Figure  43.  +4S*  Cutout  Section  Model 
Mode  1,  f  *  556  hz 


Figure  44.  +45*  Cutout  Section  Model 
Mode  2,  f  -  592  hz 


Fifure  45.  +45*  Cutout  Section  Model 
Mode  3,  f  »  711  hz 


Figure  45.  +45*  Cutout  Section  Model 
Mode  4,  f  -  753  hz 
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The  +45°  cutout  section  caused  a  loss  in  stiffness  for  all  the  modes, 
as  compared  to  full  panel  model.  This  agrees  with  the  trend  of  Goodwin’s 
work,  as  is  shown  in  Figure  49.  Both  show  a  decrease  in  the  first  natural 
frequency  and  an  increase  in  the  second,  as  compared  to  the  full  panel 
results.  Note  that  the  NISA  code  was  extremely  limited  and  unable  to 
model  this  configuration. 


Figure  49.  +45*  Cutout  Effects 


-li*  CUTOUT  SECTION  MODEL.  This  is  one  of  the  most 
interesting  models.  Its  natural  frequencies  more  closely  match  Goodwin’s 
than  any  other  cutout  (see  Table  10),  and  it  sheds  light  on  the  higher* 
frequency  mode  switching  phenomenon. 


Table  10.  -45°  Cutout  Section  Model  Results 


r, 

fj 

fy 

f4 

~~JyZ  . 

-45”  Cutout  Model 

550 

584 

698 

749 

818 

Goodwin  Exper. 

528 

538 

648 

666 

823 

%  Difference 

•t-8.55% 

+  7.72* 

+  12.5% 

-0.607* 

(*  all  frequencies  are  in  hertz) 


Its  Hrst  Hve  mode  shapes  match  those  found  by  Goodwin's 
experimental  work  as  well  as  those  found  by  the  STAGS  C-1  work  of  Cyr. 
Its  Hrst  four  mode  shapes  match  the  first  four  of  Cyr's  experimental  work, 
and  this  model's  sixth  mode  shape  matches  Cyr's  fifth  experimental  mode 
shape.  It  appears  that  Cyr  expeiienced  mode  switching  with  hi:;  fifth  mode 
but  never  noticed  it,  since  he  never  investigated  the  sixth  mode.  Rajamani 
and  Prabhakaran  noted  a  relationship  between  mode  switching  at  higher 
frequencies  and  cutout  orientation,  and  perhaps  this  model  is  more  sensitive 
to  such  behavior  due  to  its  unique  boundary  condition  [16:576].  The  order 
of  mode  shapes  presented  here  in  Figures  SO  through  55  is  thought  to  be  the 
proper  one. 

It  is  also  interesting  to  note  that  the  mode  shapes  for  the  -45°  cutout 
section  model  are  the  transposes  of  those  for  the  +45°  model,  except  that 
the  Hfth  and  sixth  mode  shape  are  switched,  as  was  previously  noted. 


Figure  54.  -45”  Cutout  Section  Model 
Mode  5,  f  s  818  hz 


Unlike  Goodwin’s  results,  this  model  shows  a  noticeable  decrease  in 
natural  frequency  for  the  first  two  mode  shapes  as  compared  to  the  full 
panel  model  as  in  shown  in  Figure  56.  However,  the  experimental 


Figure  56.  -45“  Cutout  Effects 

results  of  Cyr  also  showed  a  noticeable  decrease  which  be  attributed  to  a 
poorly  clamped  test  condition.  It  would  make  sense  for  the  natural 
frequencies  to  increase  as  they  did  for  the  +45*  model,  since  the  laminate 
is  a  quasi-isotropic  layup.  However,  the  results  for  the  -45*  orientation  are 
close  enough  to  those  of  for  the  +45*  orientation  for  the  Hrst  two  modes 
(less  than  1.5%  in  each  case)  that  this  difference  can  be  attributed  to  slight 
differences  when  defining  the  geometries  for  the  models  and  numeric  error 
in  the  solution  process. 
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REPAIRED  PANEL  MODELS 


REPAIRED  PANEL  MODEL.  This  repair  orientation  proved  the 
presence  of  mode  switching  between  the  fifth  and  sixth  mode  shapes, 
although  that  was  not  immediately  obvious.  Initially,  the  perfect  bond 
model  and  individual  rivet  model  were  solved  and  compared.  An  extra 
mode  shape  occurred  in  the  individual  rivet  model  solution  between 
Goodwin’s  fourth  and  fifth  mode  shapes.  The  mesh  was  refined  so  that  the 
elements  were  approximately  %  *  square  (see  Figure  57)  and  the  model 


rerun.  The  natural  frequencies  were  more  accurate  (see  Table  11),  but  the 
additional  mode  shape  remained.  The  fist  six  mode  shapes  for  both  models 
were  then  examined  in  close  detail. 
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Table  11.  O'*  Repair  Model  Results 


rr- 

U 

U 

. 

Rivet  Model 

R  460 

531 

702 

780 

830 

Goodwin  Exp.  (GE) 

1 

555 

663 

745 

795 

%  Difference 

1  +4.78« 

-4.32* 

+5.88* 

+4.70* 

+  4.40* 

ReHned  Rivet  Model 

II  447 

522 

689 

757 

815 

%  Diff.  vs  GE 

-5.95* 

+3.92* 

+  1.61* 

+  2.52* 

Perfect  Bond  Model 

1  537 

600 

779 

928 

954 

%  Diff.  vs  GE 

1  -^22% 

+  8.2* 

+  17* 

+24* 

+20* 

NISA 

1  490 

549 

730 

884 

N/A 

%  Diff.  vs  GE 

1  -1-11. 6« 

-1.08* 

+  10.1* 

+  18.7 

(*  all  frequencies  in  hertz) 


The  perfectly  bonded  repair  model  resulted  in  mode  shapes  that 
matched  Goodwin’s  experimental  and  FEA  results,  although  the  associated 
natural  frequencies  were  very  high.  The  individual  rivet  model,  however, 
agreed  with  Goodwin  for  the  first  four  mode  shapes,  but  showed  Goodwin’s 
fifth  mode  shape  at  the  sixth  natural  frequency.  Also,  the  fifth  mode  shape 
of  the  individual  rivet  model  matched  the  sixth  mode  shape  of  the  perfectly 
bonded  model.  This  verified  that  the  extra  mode  shape  was  the  result  of 
mode  switching.  These  individual  rivet  model  inode  shapes  are  compared 
against  the  perfect  bond  model  mode  shapes  in  Figures  57  through  69.  the 
refined  rivet  model  mode  shapes  are  shown  in  Figures  70  through  75. 
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Figure  62.  0*  Repair  Individual  Rivet  Model 
Mode  3,  f  *  702  hz 


Figure  63.  0*  Repair  Perfect  Bond  Model 
Mode  3,  f  -  779  hz 
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Figure  66.  0**  Repair  Individual  Rivet  Model 
Mode  5,  f  830  hz 


Figu-  2  67.  0*  Repair  Perfect  Bond  Model 
Mode  5,  f  «  954  hz 
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Figure  69.  0*  Repair  Perfect  Bond  Model 
Mode  6,  f  -  992  hz 
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Figure  72.  0**  Repair  Refined  Rivet  Model 
Mode  3,  f  =  689  hz 


Figure  73.  0*  Repair  Refined  Rivet  Model 
Mode  4,  f  -  757  hz 
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Figure  74.  0**  Repair  Refined  Rivet  Model 
Mode  S,  f  -  815  hz 


Figure  75.  0**  Repair  Refined  Rivet  Model 
Mode  6,  f  »  849  hz 


The  graph  in  Figure  76  show  that  individually  modeling  the  rivets 
more  closely  approximates  the  actual  behavior  of  the  panel  while  the 
assuming  a  perfect  bond  continues  the  trend  set  by  the  NISA  results  by 
being  consistently  and  increasingly  too  stiff.  This  makes  sense,  since  the 
patch  is  only  "bonded"  to  the  panel  at  the  rivets. 


Figure  76.  0"  Repair  Effects 


REEAIREP.  PANEL  MODELS.  The  natural  frequencies  of  the 

individual  rivet  model  are  a  substantial  improvement  over  the  perfect  bond 

♦ 

model  and  a  good  approximation  of  the  actual  conditions,  as  shown  in 
Table  12.  Both  the  individual  rivet  model  and  the  perfect  bond  model  show 
mode  switching  between  modes  four  and  five  with  respect  to  Goodwin’s 
experimental  work.  They  agree  with  the  mode  shapes  as  predicted  by 


Table  12.  90’’  Repair  Model  Results 


fr 

fr 

f4 

f5 

1  Rivet  Model 

475 

506 

710 

759 

828 

I  Goodwin  Ezp.  (G£) 

509 

523 

706 

758 

796 

1  %  Difference 

•6.67% 

-3.25% 

+0.566% 

+0.132% 

+  4.02%, 

1  Perf.  Bond  Model 

592 

616 

817 

955 

978 

1  *  Diff.  vs  GE 

-f-16% 

+  13% 

.  +16% 

+  26% 

+  23% 

NISA 

500 

535 

749 

800 

N/A 

%  Diff.  vs  GE 

-1.77% 

+2.29% 

+6.09% 

+5.54% 

1 

1 

1 

(*  all  frequencies  are  in  hertz) 


Goodwin’s  NISA  work,  and  they  are  consistent  with  the  full  panel  mode 
shapes.  These  circumstances  show  that  Goodwin  experienced  mode 
switching  between  his  fourth  and  fifth  mode  shapes  in  his  experimental 
work.  The  mode  shapes  shown  here  are  therefore  felt  to  be  ordered  in  the 
proper  arrangement.  The  mode  shapes  are  presented  in  Figures  77  through 


86. 


Fisure  77.  90“  Repair  Individual  Rivet  Model 
Mode  1,  f  =  475  hz 


Figure  78 


90“  Repair  Perfect  Bond  Model 
Mode  1.  f  -  592  hz 


Figure  85. 


90'’  Repair  Individual  Rivet  Model 
Mode  S,  f  »  828  hz 


Figure  86.  90*  Repair  Perfect  Bond  Model 
Mode  5,  f  *  978  hz 
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Figure  87  allows  a  graphic  ccmparison  of  the  frequency  results.  The 
most  striking  feature  is  that  the  first  tv/o  nstur'’.!  fnequencies  of  the 
individual  rivet  model  were  lower  than  the*  firs-  two  natural  frequencies 
reported  by  Goodwin  for  his  experimental  worl'.  This  is  unusun;  because 
the  models  using  QUAD4  elements  have  been  unnaturally  aiiff  before  this. 


Figure  87.  90*  Repair  Effects 


Since  this  90*  orientation  causes  the  patch  to  cross  more  node  lines  than 
the  0*  orientation,  tne  presence  of  the  CONM2  elements  at  the  rivet 
locations  could  be  responsible  for  this  drop  in  the  first  two  natural 
frequencies.  These  elements  add  point  masses  to  these  locations  without 
adding  any  stiffness.  The  90*  orientation  allows  for  more  node  lines  to  be 
affected.  At  the  higher  frequencies,  the  node  li.ies  are  more  sparsely 
scattered  across  the  patch  region,  which  would  lessen  this  effect. 
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REPAIRED  PANEL  MODEL.  These  models  behaved  very 
much  as  expected,  and  they  are  very  similar  to  the  next  set,  as  would  be 
expected  from  symmetric  repairs  to  a  quasi-isotropic  laminate.  The 
frequency  results  are  shown  in  Table  13. 


Table  13.  +45”  Repair  Model  Results 


fr 

f4 

fs 

Rivet  Model 

510 

737 

799 

854 

Goodwin  Exp.  (GE)  | 

506 

525 

721 

789 

820 

%  Difference 

-3.56ft 

-2.86ft 

+2.22ft 

+  1.27* 

+4.15* 

Perfect  Bond  Model 

559 

583 

783 

886 

917 

%  tuff,  va  GE 

+  10.5ft 

+  11.0ft 

8.60ft 

+  13.3* 

+  11.8*  1 

(*  all  frequeaciea  are  ia  hertz) 


The  mode  shapes  themselves  are  also  well-behaved,  except  for  a 
mode  switch  between  the  fourth  and  fifth  modes  of  both  of  these  models 
with  respect  to  Goodwin's  experimental  mode  shapes.  Since  the  mode 
shapes  shown  in  Figures  88  through  97  match  those  of  the  undamaged 
panel,  they  are  believed  to  be  arranged  in  the  proper  order. 
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Figure  90.  +45**  Repair  Individual  Rivet  Model 
Mode  2,  f  «  SIO  hz 


Figure  91.  +45*  Repair  Perfect  Bond  Model 
Mode  2,  f  -  583  hz 
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Figure  94.  +45”  Repair  Individual  Rivet  Model 
Mode  4,  f  -  799  hz 


Figure  95.  +45*  Repair  Perfect  Bond  Model 
Mode  4,  f  -  886  hz 
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Figure  96.  +45®  Repair  Individual  Rivet  Model 
Mode  5,  f  854  hz 


Figure  97.  +45*  Repair  Perfect  Bond  Model 
Mode  5,  f  -  917  hz 


The  graph  in  Figure  98  show  that  for  the  first  two  natural 
frequencies  the  individual  rivet  model  are  lower  than  those  reported  by 
Goodwin  for  his  experimental  work. 


UadeNimfaer 

SBC**- SSrmhmw  S*«f.  mim 


Figure  98.  +45*  Repair  Effects 

This  is  the  same  behavior  exhibited  by  the  90^  repaired  panel  model, 
though  to  a  lesser  extent.  Since  this  patch  is  oriented  halfway  between  the 
0**  and  90**  patch,  the  behavior  of  this  model  should  reflect  aspects  of  both. 

-12**  REPAIRED  PANEL  MODEL.  This  results  from  this  model  are 
very  similar  to  the  preceding  one,  as  previously  stated.  .  The  frequency 
results  are  shown  in  Table  14. 
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Table  14.  -45*"  Repair  Model  Results 


u 

f* 

f5 

- 

Rivet  Model 

476 

513 

. 

702 

775 

842 

Goodwin  Exper.  (GE) 

505 

516 

717 

750 

806 

%  Difference 

-5.74* 

-0.581  * 

-2.09* 

+  3.33* 

Perfect  Bond  Model 

608 

642 

831 

943 

948 

%  Diff.  V8  GE 

+  20* 

+24* 

+  16* 

+26* 

+  18* 

(*  til  frequeocies  tre  in  hertz) 


The  mode  shapes  for  the  individual  rivet  model  agree  with  those  from 
Goodwin’s  experimental  work,  but  the  fourth  and  fifth  mode  shapes  of  the 
perfect  bond  model  are  indistinct  and  cannot  be  accurately  identiHed.  The 
mode  shapes  are  shown  in  Figures  99  through  108. 


93 


Figure  99.  -45**  Repair  Individual  Rivet  Model 
Mode  1,  f  =  476  hz 


Figure  100.  -45  *  Repair  Perfect  Bond  Model 
Mode  1,  f  •  608  hz 
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Figure  IGl.  -45°  Repair  Individual  Rivet  Model 
Mode  2,  f  »  513  hz 


Figure  102.  -45°  Repair  Perfect  Bond  Model 
Mode  2,  f  »  642  hz 
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Figure  103.  -45°  Repair  Individual  Rivet  Model 
Mode  3,  f  =  702  hz 


Figure  104.  -45°  Repair  Perfect  Bond  Model 
Mode  3,  f  «  831  hz 
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Figure  107.  >45*  Repair  Individual  Rivet  Model 
Mode  5,  f  =  S42  hz 


Figure  108.  -45*  Repair  Perfect  Bond  Model 
Mode  5,  f  -  948  hz 
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The  graphic  comparison  shown  in  Figure  109  is  very  similar  to 
Figure  98,  as  expected. 


Figure  109.  -45®  Repair  Effects 
The  individual  rivet  model  follows  the  trend  of  the  experimental 
work  by  Goodwin,  whereas  the  perfect  bond  model  is  consistently  high. 
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V  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  has  shown  that  a  curved  rectangular  composite  panel 
modeled  with  flat  MSC/NASTRAN  QUAD4  elements  can  be  successfully 
modifed  to  show  the  effect  of  a  cutout  region  and  subsequent  repair  on  the 
dynamic  behavior  of  the  panel.  More  specifically,  using  ELAS2  and 
CONM2  elements  to  model  the  individual  rivets  in  conjunction  with  GAP 
elements  to  model  the  separation  of  the  midsurfaces  of  the  panel  and  the 
repair  allowed  a  much  more  accurate  prediction  of  the  first  Hve  natural 
frequencies  of  the  repaired  panel  than  did  a  model  assuming  a  perfect  bond 
between  the  patch  and  the  panel.  The  individual  rivet  model  also  was 
better  able  than  the  perfect  bond  model  to  predict  accurately  the  first  five 
mode  shapes  of  the  repaired  panel,  with  the  exception  of  the  occasional 
switch  between  the  fourth  and  fifth  mode  shapes.  It  was  able  to  show  the 
effect  of  the  rivets  on  the  mode  shapes,  whereas  the  perfect  bond  model 
was  not. 

The  occurance  of  mode  switching  in  these  models  is  not  believed  to 
be  important,  though.  The  mode  switching  only  occurs  between  the  higher 
mode  shapes,  which  are  of  less  relevance  than  the  first.  Also,  since  mode 
switching  is  experienced  in  experimental  models  as  well  and  analytical 
ones,  it  is  likely  to  be  encountered  outside  the  laboratory  and  should  not  be 
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ignored. 

For  further  research  in  this  area,  the  static  failure  modes  of  a 
repaired  panel  should  be  investigated,  both  experimentally  and  analytically. 
Along  a  similar  vein,  the  effect  of  the  riveted  repair  on  the  buckling 
behavior  of  the  repaired  panel  should  be  investigated. 

It  is  unlikely  that  the  future  of  ABDR  will  rest  with  riveted  repairs. 
The  advances  being  made  with  low  temperature,  fast-curing  adhesives  make 
it  seem  probable  that  they  will  see  more  use  in  field  repairs.  Therefore, 
the  dynamic  behavior  of  such  a  repair  should  be  investigated,  as  well  as  its 
failure  modes  and  the  its  effect  on  the  buckling  behavior  of  the  panel. 
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